####Boxplot of contributions
setwd('/media/Iomega_HDD/boulot/Article/kd/last_23_01_2012/data/')
load("dataset.Rdata")
table <- cbind(signif(t5$value,7), signif(t5$gsm_chl,7), signif(t5$gsm_cdm,7),
                   signif(t5$gsm_bbp,7),t5$lambda, t5$convert_lambda)

para<-read.table('../script_admb/model.std',header=T, fill=T)
para <- para[,1:4]

s0<-para[9,3]
aphstar<-para[27:34,3]
a2<-para[18:25,3]
eta<-para[26,3]


dl<-matrix(0,nrow=1,ncol=3)
for (i in 1:nwave){
  kds<-table[which(table[,6]==i),1]		
	chl_gsm<-table[which(table[,6]==i),2]	
	acdm_gsm<-table[which(table[,6]==i),3]
	bbp_gsm<-table[which(table[,6]==i),4]	
	S<-s0
	CC<-acdm_gsm
	BBP<-bbp_gsm
	
	Kdmod<-kw[i]+CC*exp(-S*(wave[i]-443))+aphstar[i]*chl_gsm^a2[i]+BBP*(wave[i]/443)^(eta)
	titi<-cbind(kds,Kdmod, wave[i])
	dl<-rbind(dl,titi)
}
dl<-dl[2:length(dl[,1]),]

dl2<-matrix(0,nrow=1,ncol=5)
for (i in 1:nwave){
	a<-dl[which(dl[,3]==wave[i]),2]
	kwu<-kw[i]
	chl_gsm<-table[which(table[,6]==i),2]	
	acdm_gsm<-table[which(table[,6]==i),3]
	bbp_gsm<-table[which(table[,6]==i),4]
	S<-s0
	CC<-acdm_gsm
	BBP<-bbp_gsm
  adm<-CC*exp(-S*(wave[i]-443))
  ch<-aphstar[i]*chl_gsm^a2[i]
  bp<-BBP*(wave[i]/443)^(eta)
  #a<-kw[i]+adm+ch+bp
  tutu<-cbind(kwu/a,adm/a,ch/a,bp/a,rep(wave[i],length(a)))
  dl2<-rbind(dl2,tutu)
}
dl2<-dl2[2:length(dl2[,1]),]

postscript('../Figures/Figure6.ps', width=14,height=14,paper="special", horizontal=F)
opar<-par()
#liml<-c(1,1,1,0.3)
par(mfrow=c(2,2))
par(mar=c(5,5,3,1))
titre<-c(expression(paste(K[w],"/",K[d], sep="")),expression(paste(a[CDM],"/",K[d])),expression(paste(a[chl],"/",K[d])),expression(paste(b[bp],"/",K[d])))
texti<-c("(a)","(b)","(c)","(d)")
for (i in 1:4){
boxplot(dl2[,i]*100~dl2[,5], xlab=expression(paste(lambda, "(nm)")), ylab="Contribution to Kd (%)", main=titre[i], cex.main=2, cex.axis=2, cex.lab=2 , ylim=c(0,100))
text(1,95, labels=texti[i], cex=2)
}
par(opar)
dev.off()
